/*
 *  hyperelasticInvariant.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 1/04/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */
#include "hyperelasticinvariants.h"

HyperelasticInvariants::HyperelasticInvariants()
:	Material(),
 	Wvol(),
	W1(),
	W2(),
	W3(),
	dWvol(),
	Wiso1(),
	Wiso2(),
	W11(),
	W12(),
	W13(),
	W22(),
	W23(),
	W33(),
	ddWvol(),
	Wiso11(),
	Wiso12(),
	Wiso22(),
	delta1(),
	delta2(),
	delta3(),
	delta4(),
	Wiso111(),
	Wiso112(),
	Wiso122(),
	Wiso222(),
	Wiso221(),
	Wiso211(),
	kappa11(),
	kappa12(),
	kappa21(),
	kappa22(),
	kappa31(),
	kappa32(),
	kappa41(),
	kappa42()
{}

// ***

HyperelasticInvariants::HyperelasticInvariants(const string& na, const string & t)
:	Material(na,t),
 	Wvol(0.0),
	W1(0.0),
	W2(0.0),
	W3(0.0),
	dWvol(0.0),
	Wiso1(0.0),
	Wiso2(0.0),
	W11(0.0),
	W12(0.0),
	W13(0.0),
	W22(0.0),
	W23(0.0),
	W33(0.0),
	ddWvol(0.0),
	Wiso11(0.0),
	Wiso12(0.0),
	Wiso22(0.0),
	delta1(0.0),
	delta2(0.0),
	delta3(0.0),
	delta4(0.0),
	Wiso111(0.0),
	Wiso112(0.0),
	Wiso122(0.0),
	Wiso222(0.0),
	Wiso221(0.0),
	Wiso211(0.0),
	kappa11(0.0),
	kappa12(0.0),
	kappa21(0.0),
	kappa22(0.0),
	kappa31(0.0),
	kappa32(0.0),
	kappa41(0.0),
	kappa42(0.0)
{}

// ***

void HyperelasticInvariants::calcInvariants( const VectorXd & C)
{
	InverseVoigtFormat(CSym,C);
	
	I = CSym.trace();
	
	II = 0.5 * (I * I - (CSym * CSym).trace());
	
	III = CSym.determinant();
	
	J = sqrt(III);
	
	if (J<=0) cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
	
	if(!hypertype.compare("decoupled"))
	{
		Jm =  pow(J, - 2.0/3.0);
		Im = Jm * I;
		IIm = Jm * Jm * II;

	}
	
	//cout << "Ci-1= " << endl << CiTensor << endl;

//  cout << "I-1= " << I << endl;
//  cout << "II-1= " << II << endl;
//  cout << "III-1= " << III << endl;
//  cout << "Jm-1= " << Jm << endl;

}

// ***

void HyperelasticInvariants::calcCTensor (const VectorXd & C)
{
    InverseVoigtFormat(CTensor,C);
}

// ***

void HyperelasticInvariants::calcCiTensor( const VectorXd & C)
{
	//InverseVoigtFormat(CTensor,C);
    InverseVoigtFormat(CSym,C);

 	CiSym = CSym.inverse();

	CiTensor(1,1) = CiSym(0,0);
	CiTensor(2,2) = CiSym(1,1);
	CiTensor(3,3) = CiSym(2,2);
	CiTensor(2,1) = CiSym(0,1);
	CiTensor(3,1) = CiSym(0,2);
	CiTensor(3,2) = CiSym(1,2);
}

// ***

void HyperelasticInvariants::calcP(const VectorXd & C)
{
	InverseVoigtFormat(CTensor,C);
    calcCiTensor(C);
    P = Id4 - 1.0 / 3.0 * dyadic(CiTensor, CTensor);

}

// ***

void HyperelasticInvariants::calcPt(const VectorXd & C)
{
	InverseVoigtFormat(CTensor,C);
    calcCiTensor(C);
    Pt = Id4 - 1.0 / 3.0 * dyadic(CTensor, CiTensor);
}

// ***

void HyperelasticInvariants::calctildeP(const VectorXd & C)
{
    calcCiTensor(C);
    tildeP = selfOdot(CiTensor) - 1.0/3.0 * selfDyadic(CiTensor);
}

// ***

void HyperelasticInvariants::calcW( const VectorXd & C )
{
	calcInvariants(C);
	strainEnergy();
}

// ***

void HyperelasticInvariants::calcS( const VectorXd & C)
{
	calcInvariants(C);	
	DstrainEnergy();
   
	if(!hypertype.compare("coupled"))
	{
		// S = 2.0 * [ (W1 + I * W2) * 1 - W2 * C + III * W3 * C^{-1}]
		
		SMat = 2.0 * ((W1 + I * W2) * Matrix3d::Identity() - W2 * CSym + III * W3 * CSym.inverse());
		VoigtFormat(SMat,S);
	}
	else
	{
		calcSiso(C);
		calcSvol(C);
		S = Svol + Siso;
	}
}

//***

void HyperelasticInvariants::calcSvol( const VectorXd & C)
{
	calcInvariants(C);
	DstrainEnergy();
	
	SvolMat = J * dWvol * CSym.inverse();
	VoigtFormat(SvolMat, Svol);
}

//***

void HyperelasticInvariants::calcSiso( const VectorXd & C)
{
	calcInvariants(C);
	DstrainEnergy();
    mSisoTens = 2.0 * ((Wiso1 + Im * Wiso2) * Id - Wiso2 * Jm * CTensor);

    calcP(C);
    //cout << "P= " << endl << P << endl;
    SisoTens = Jm * doubleContraction(P, mSisoTens);
    
	VoigtFormat(SisoTens, Siso);
}


//***

void HyperelasticInvariants::calcCt( const VectorXd & C)
{
	if(type.compare("svk"))
	{
		//int i,j,k,l;
		calcInvariants(C);
		DstrainEnergy();
		DDstrainEnergy();
		
		if(!hypertype.compare("coupled"))
		{  
			delta1 = 4.0 * (W2 + W11 + 2.0 * I * W12 + I * I * W22);
			delta2 = 4.0 * (W12 + I * W22);
			delta3 = 4.0 * (III * W13 + I * III * W23);
			delta4 = 4.0 * (III * W3 + III * III * W33);
            
            calcCiTensor(C);
			//cout << "Ci-2= " << endl << CiTensor << endl;
			IIc = selfOdot(CiTensor); // In this function is computed both CTensor, CiTensor and IIc
      
			CtTens = delta1 * selfDyadic(Id) - 4.0 * W2 * Id4 - 4.0 * III * W3 * IIc - delta2 * ( dyadic(Id,CTensor) + dyadic(CTensor,Id))
			    + 4.0 * W22 * selfDyadic(CTensor) + delta3 * (dyadic(Id,CiTensor) + dyadic(CiTensor,Id)) 
			    - 4.0 * III *  W23 * (dyadic(CTensor,CiTensor) + dyadic(CiTensor,CTensor)) + delta4 * selfDyadic(CiTensor);
			
			VoigtFormat(CtTens, Ct);
		}
		else // Hyperelastic Material defined in decoupled form.
		{  
			delta1 = 4.0 * (Wiso11 + 2.0 * Im * Wiso12 + Wiso2 + Im * Im * Wiso22);
			delta2 = - 4.0 * (Wiso12 + Im * Wiso22);
			delta3 = 4.0 * Wiso22;
			delta4 = - 4.0 * Wiso2;
			
//			CiSym = CSym.inverse();
            
            calcCiTensor(C);

			IIc = selfOdot(CiTensor);
			
			//First, the Volumetric Part of Ct is filled: Ct_vol = J(p + pdp/dJ)Ci * Ci - 2Jp dC^-1/dC
			
			CtTens = J * (dWvol + J * ddWvol) * selfDyadic(CiTensor) - 2.0 * J * dWvol * IIc;
			
			VoigtFormat(CtTens,Ctvol);
			
// 			cout << "Ctvol= " << endl << Ctvol << endl << endl;
			
			//Holzapfel's book pag. 255: expression for Ctm_Iso tensor
			
			//First, it is needed the modified Constitutive Tensor, 
			
			P = Id4 - 1.0 / 3.0 * dyadic(CiTensor, CTensor); //getP(C) but we already know CTensor
			Pt = Id4 - 1.0 / 3.0 * dyadic(CTensor, CiTensor); //getPt(C)
			
			mCtisoTens = Jm * Jm * (delta1 * selfDyadic(Id) + delta2 * (dyadic(Id, Jm * CTensor) + dyadic(Jm * CTensor,Id)) + delta3 * selfDyadic(Jm * Jm * CTensor)
			        + delta4 * Id4);
			
			//Computing first term 
			CtisoTens = doubleContraction(P,(doubleContraction(mCtisoTens, Pt)));
			
			//Computing second term of Ctm_Iso
			CtisoTens += 2.0 / 3.0 * Jm * doubleContraction(mSisoTens,CTensor) * (IIc - selfDyadic(CiTensor) / 3.0);
			
			//Computing third term:
			CtisoTens -= 2.0 / 3.0 * (dyadic(CiTensor,SisoTens) + dyadic(SisoTens,CiTensor));
			
			VoigtFormat(CtisoTens,Ctiso);
			
			//To conclude, the sum of the both parts computed above is performed in order to obtain the total elasticity tensor:
			
			CtTens += CtisoTens;
			
			VoigtFormat(CtTens, Ct);
		}
	}
    
   // cout << "Ct_iso= " << endl << Ctiso << endl << endl;
}

//***
